bzoj 4816 [Sdoi2017]数字表格.md

题目链接

bzoj4816: [Sdoi2017]数字表格

题解

满满的反演的套路的味道

$$Ans= \prod_{i=1}^{n} \prod_{j=1}^{m} f[gcd(i,j)]$$

常规操作枚举约数

$$\prod_{d=1}^n\prod_{i=1}^n\prod_{j=1}^m gcd(i,j)==d\ ?\ f[gcd(i,j)]$$

上面的式子可以化为

$$\prod_{d=1}^{n}f[d]^{\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]}$$

对于右上角的式子

$$\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]$$

很眼熟…推下式子

$$=\sum_{i=1}^{? \dfrac{n}{d}?}\sum_{j=1}^{?\dfrac{m}{d}?}[gcd(i,j)==1]$$ $$=\sum_{i=1}^{min(\dfrac{n}{d},\dfrac{m}{d})} \mu(i)\sum_{i|p,p\leq n/d }\sum_{i|j,j\leq m/d}1$$ $$=\sum_{i=1}^{n/d}\mu(i)\left \lfloor \frac{n}{id}\right \rfloor \left \lfloor\frac{m}{id} \right \rfloor$$

令$p=id$那么原式变成了

$$\prod_{p=1}^{min(n,m)} \prod_{d|p} f[d]^{\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor \mu(\frac{p}{d})} = \prod_{p=1}^{min(n,m)} (\prod_{d|p} f[d]^{\mu(\frac{p}{d})})^{\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor}$$

对于每个$p$,令$F[p]=\prod_{d|p} f[d]^{\mu(\frac{p}{d})}$,$F[p]$的值是固定的，用筛法求出$F[p]$，做前缀积，对与$\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor$数论分块一下
复杂度为$O((n+T\sqrt{n})log(1e9+7))$

代码

#include<cstdio>
#include<algorithm>
const int maxn = 1000007;
const int mod = 1e9+7;
int f[maxn+7],invf[maxn+7],F[maxn+7];
bool isprime[maxn+7];int mu[maxn+7],prime[maxn+7],num;

int pow(int a,int p) {
    int ret=1;
    for(;p;p>>=1,a=(1LL*a*a)%mod){
        if(p&1) ret=(1LL*ret*a)%mod;
    }
    return ret;
}

void init() {
    f[1]=1;
    for(int i=2;i<maxn;++i) f[i]=(f[i-1]+f[i-2])%mod;
    for(int i=1;i<maxn;++i) {
        invf[i]=pow(f[i],mod-2)%mod;
        //printf("%d %d\n",invf[i],f[i]);
    }
    isprime[1]=1;mu[1]=1;F[0]=F[1]=1;
    //printf("%d\n",mu[1]);
    for(int i=2;i<=maxn-7;++i) {
        F[i]=1;
        if(!isprime[i]) prime[++num]=i,mu[i]=-1;
        for(int j=1;j<=num&&prime[j]*i<=maxn-7;++j) {
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=maxn-7;++i) {
        if(mu[i]==0)continue;
        for(int j=i;j<=maxn-7;j+=i) {
            if(mu[i]==1) F[j]=(1ll*F[j]*f[j/i])%mod;
            if(mu[i]==-1) F[j]=(1ll*F[j]*invf[j/i])%mod;
        }
    }
    for(int i=1;i<=maxn-7;++i) F[i]=(1ll*F[i]*F[i-1])%mod;//,printf("%d\n",F[i]);
}
int query(int a,int b) {
    int n=std::min(a,b);long long ans=1;
    for(int nxt,i=1;i<=n;i=nxt+1) {
        nxt=std::min(a/(a/i),b/(b/i));
        long long tmp=1ll*F[nxt]*pow(F[i-1],mod-2)%mod;
        ans=(1ll*ans*pow(tmp,1ll*(a/i)*(b/i)%(mod-1)))%mod;
    }
    return (ans+mod)%mod;
}

int main() {
    //freopen("001.out","w",stdout);
    init();	
    int T;scanf("%d",&T);
    for(int a,b;T--;) {
        scanf("%d%d",&a,&b);
        printf("%d\n",query(a,b));
    }
    return 0;
}